

function vesselID = manualSegmentArteries(avios,traj,numarteries)
% allows you to manually segment arteries for the high res widefield
% tracking

vesselID = nan(size(traj(1).vesselCent,1),1);
SE = strel("disk",15);


plotfig = figure;
imagesc(avios);colormap('gray');axis equal;
hold on
plot(traj(1).vesselCent(:,1),traj(1).vesselCent(:,2),'r.');

if nargin < 3 || isempty(numarteries)
    numarteries = input('how many arteries? ');
end

n = 0;
while n < numarteries
    selectionLine = drawpolyline;
    sent = input('keep line? [1/0]: ');
    if sent ~= 1
        continue;
    end
    n = n + 1;
    startpoint = false(size(avios));startpoint(round(selectionLine.Position(1,2)),round(selectionLine.Position(1,1))) = true;
    BW = createMask(selectionLine);
    BW = imdilate(BW,SE);
    D = bwdistgeodesic(BW,startpoint);
    temp = round(traj(1).vesselCent);
    temp = sub2ind(size(avios),temp(:,2),temp(:,1));
    temp(~isnan(temp)) = D(temp(~isnan(temp)));
    %temp = D(temp);
    
    [distanceAlongVasc,pointidx] = sort(temp,'ascend');
    pointidx = pointidx(~isnan(distanceAlongVasc));    
    plot(traj(1).vesselCent(pointidx,1),traj(1).vesselCent(pointidx,2),'r.-')
    vesselID(pointidx) = n;
end

figure(plotfig)
hold off
imagesc(avios);colormap('gray');axis equal;
hold on
clrs = jet(numarteries);
x = traj(1).vesselCent(~isnan(vesselID),1);
y = traj(1).vesselCent(~isnan(vesselID),2);
z = clrs(vesselID(~isnan(vesselID)),:);
scatter(x,y,15,z,'filled');
hold off


